arXiv: 1506.01193vl [math.NA] 3 Jun 2015 


Locally Supported Wavelets for the Separation of Spherical 
Vector Fields with Respect to their Sources 


C. Gerhards 

Abstract We provide a space domain oriented separation of magnetic fields into parts generated 
by sources in the exterior and sources in the interior of a given sphere. The separation itself 
is well-known in geomagnetic modeling, usually in terms of a spherical harmonic analysis or a 
wavelet analysis that is spherical harmonic based. In contrast to these frequency oriented meth¬ 
ods, we use a more spatially oriented approach in this paper. We derive integral representations 
with explicitly known convolution kernels. Regularizing these singular kernels allows a multiscale 
representation of the internal and external contributions to the magnetic held with locally sup¬ 
ported wavelets. This representation is applied to a set of CHAMP data for crustal held modeling. 
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1 Introduction 

The Earth’s magnetic held is a complex structure consisting of various contributions, such as the 
dominating core held, the crustal held, and effects from iono- and magnetospheric processes. A 
major task in understanding the geomagnetic held is the separation of these contributions. An 
overview on different approaches to this is given, e.g., in |24| . A hrst step is the mathematical sep¬ 
aration of magnetic held measurements taken at satellite altitude into contributions from sources 
in the exterior of the orbit and contributions from sources in the interior. Generally, we assume 
the magnetic held b to be governed by the pre-Maxwell equations 

V A 6 = noj, 

V • 6 = 0, 

with j describing the source current density and /tq the vacuum permeability (A denotes the vector 
product). If no source currents j are present, one has b = VC/, for some harmonic potential U, 
and the typical approach to modeling the magnetic held is the so-called Gauss representation of 
the corresponding potential in terms of scalar spherical harmonics Yn^k (see, e.g., [3] and [25]). 
Generally, however, satellite data is collected in a source region of the magnetic held. Then the 
Mie decomposition allows a decomposition of the magnetic held into a poloidal part pb and a 
toroidal part qb- The toroidal part describes the magnetic held due to poloidal current densities 
Pj , while the poloidal part can be split into a part p^^^ that is due toroidal sources in the exterior 
of the satellite’s orbit, and a part that is due to toroidal sources in the interior. A more 
detailed description can be found, e.g., in [2] and |3]. In this setting, the quantities pf^* and 
qb can be expanded in a system of vector spherical harmonics \ respectively (a 

system that actually originates in quantum mechanics; see, e.g., m- ’ 

However, due to the global nature of scalar and vector spherical harmonics, they are not the best 
choice for modeling strongly localized structures, such as the Earth’s crustal held, or modeling 
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from only locally available data. Several multiscale approaches with spatially better localizing 
kernels have been developed to improve this drawback, e.g., in [5], m for potential fields, in 
[55], [53] for the above described separation with respect to the sources, and in |3], [5D] for a 
representation of ionospheric magnetic fields and current densities. A comprehensive introduction 
of kernel functions for such methods can also be found in [T5] . 

It is the aim of this paper to transfer the multiscale approach described in [55] , [53] , which is based 
on a construction of scaling and wavelet kernels in frequency domain (i.e., based on an adequate 
superposition of the vector spherical harmonics y^^), to a setting where the scaling and wavelet 
kernels are constructed entirely in space domain. For that purpose, the vector spherical harmonics 
fc> * = 1) 3, are described by operators z = 1, 2, 3. A decomposition of the magnetic field 

in terms of these operators, in combination with the spherical Helmholtz decomposition, allows 
an integral expression of the quantities Qb- Motivated by [T3], a regularization of the 

convolution kernels appearing in this integral expression provides a multiscale representation with 
wavelets that are locally supported in space. The multiscale representation is described in detail 
in Section There, we also apply the derived algorithm to a set of real CHAMP satellite data. 
The preparatory construction of the regularized kernels and a decomposition with respect to the 
operators 0 ^*^ is described in Sections]^ and Section[^provides fundamental aspects on Legendre 
polynomials and scalar and vector spherical harmonics. 


2 Preliminaries 


By Pn : [—1,1] —t n S No, we denote the set of Legendre polynomials of degree n, by 

Yn^k : H —)■ K, n e No, k = 1,... ,2n + 1, an orthonormal set of spherical harmonics of degree n 
and order k {flu = {cc G K^| |a;| = i?} denotes the sphere of radius i? > 0 and H = Hi the unit 
sphere). The fundamental connection between these two function systems is the so-called addition 
theorem, 

2n+l 2 I 1 

Y.^r.,kmnAri) = ^^Pn{£.-ri). ^ G 

k=l 


This allows us to expand zonal kernels (i.e., functions F : H x H —>■ K that satisfy F(^, rf) = G{^-r]), 
^,77 G H, for an adequate function G : [—1,1] —t K) in terms of Legendre polynomials. Known 
closed representations for certain series of Legendre polynomials can then be used to derive closed 
representations for some zonal kernels appearing in this paper. One of these series is the generating 
series for the Legendre polynomials. 


Y^h^Pnit) 

n—0 


1 

Vl + 


t G [—1, 1], h G (—1, 1). 


From this, one can derive various further representations that are, e.g., listed in [16] . Of importance 
to us are the following ones. 


Lemma 2.1. For t G (—1,1), we have 


Y -Pn{t) = In 

^ 71 


Y 


n + l 




Pnit) = In 1-f 


1-P 

VT^t 


+ In (2) , 


- 1 . 


Furthermore, the generating series for the Legendre polynomials yields an expansion of the single 
layer kernel. This is of interest since it allows an integral definition of the single layer operator 
and a definition in terms of pseudodifferential operators. 
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Lemma 2.2. Let x,y with jccj < \y\. Then 


^ ao 

^ =^y 
\^-y\ l2^l ^ 


\y\ 


3^ ^ ^ 

|a:| I 2 /I 


The set of the previously mentioned spherical harmonics yields a complete orthonormal system in 
L^(n) = {F : n —>■ Kj \F{ri)\'^du!{r]) < 00 }. The modeling of magnetic fields, however, is in first 
place a vectorial problem. For that purpose, we introduce two different complete sets of vector 
spherical harmonics. The first set requires the operators 

of^FiO = mo, (2.1) 

ofF(0 = Cell, (2.2) 

ofFiO = L^FiO, ( 2 . 3 ) 

for F : n —)■ K a sufficiently smooth scalar function, V* the surface gradient (i.e., the tangential 
part of the gradient V; more precisely, for x = G with r = |a;|, ^ = |^), 

and L* the surface curl gradient (acting as V| in a point ^ G Lt). A complete orthonormal 

system in P{Lt) = {/ : n —>■ lf(y)l'^duj(y) < 00 } is then given via 


yn,k = iTn^) i = 1,2,3, n e No,, fc = 1,... ,2n +1, 


(2.4) 


( 2 ) 

where 0^ is an abbreviation for 0i = 0 and 0^ = 1, i = 2,3, and denotes the normalization 
constants fin'^ = 1 and fXn'^ = n(n + !),'« = 2, 3. Concerning the notation, upper case letters, such 
as F, generally denote scalar valued functions, lower case letters, such as /, denote 

vector valued functions, and bold face letters denote tensor valued functions. The same notation 
holds for the function spaces C'^^^(n), of fc-times continuously differentiable functions and 

the spaces L'^{Ll), P{Ll) of square integrable functions. 

The second set of vector spherical harmonics requires the modified operators 


oW = + 


( 2 . 5 ) 

5 ( 2 ) = 


( 2 . 6 ) 

II 


( 2 . 7 ) 

0 =(-A- + i)*- 


( 2 . 8 ) 


where 


By A* we denote the Beltrami operator V* • V*. The operator D is treated in more detail in 
Subsection 3.2 A second complete orthonormal system in P{n) is then given via 


Vni = z = 1, 2,3, n e No„ fc = 1,..., 2n + 1, 


( 2 . 9 ) 


where fln'^ denotes the normalization constants = [n + l)(2n + 1), = n{2n + 1) and 

/qN 

fin = n{n + 1). The advantage of this basis system is its connection to the inner and outer 
harmonics, i.e., the functions H^Kx) = (|f|), x G = {x G |a;| < R}, and 

^Tki^) — Ti{]x^'^^^Yn,k{-y), X G flff* = {x G |a;| > i?}, which yield solutions to the inner 
and outer Dirichlet boundary value problem, respectively (i.e., boundary values Fl™l = FL'^^ = 
Yn,k on TIr and = 0 in AiF“^* = 0 in ‘). We have 

T=\xlx = rS,G^\ ( 2 . 10 ) 


( 1 ) 


-,( 2 ) 


-VxHl%{x) = 


/ d \ ^+2 _ 

(7; (An^)®yi^fc(C), r =\x\,x = riGLy. ( 2 . 11 ) 
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For a more comprehensive introduction of the function systems mentioned in this section, the 
reader is referred to, e.g., [H] and the references therein. The special importance of the last set 
of vector spherical harmonics in geomagnetic modeling is well emphasized, e.g., in m, m and 
[23| . In this paper, however, they are only to be understood as a motivation for the Helmholtz 
decomposition and a modified decomposition with respect to 0 *-*^ Our main goal is to actually 
avoid spherical harmonic representations. 


3 Regularized Kernels 

Green’s function for the Beltrami operator and the single layer kernel are especially useful when 
working with differential equations involving the operators V*, L*, A* and D. We briefly re¬ 
capitulate some of the properties of these functions and the corresponding operators before we 
introduce a regularization for both kernels separately and for their combination. To achieve inte¬ 
gral representations for the scalars of the classical Helmholtz decomposition, it is actually sufficient 
to only have Green’s function. The single layer kernel becomes necessary when we introduce a 
decomposition that pays tribute to interior and exterior sources. 

3.1 Green’s Function 

By Green’s function with respect to the Beltrami operator we denote the uniquely defined function 
G(A*; •) : 1,1) —>■ M satisfying the properties 

(i) r] I—G(A*; ^ • ry) is twice continuously differentiable on the set {ry € 1 — ^ • ry > 0}, and 

A;G(A*;C-ry) = -^, 1 - ^ ry > 0, 

for any fixed ^ € H, 

(ii) for any fixed ^ S H, the function 

ry G{A*-f ' v) - ^ ln(l 

is continuously differentiable on O, 

(hi) for any fixed ^ € H, 

^ ■ p)duj{ri) = 0. 

One can verify the following explicit representation, 

G(A*;^,y) = ^ ln(l - ^ • ry) + ^(1 - ln(2)), 1 - ^ • ry > 0. (3.1) 

The bilinear series expansion reads 

oo 2n+l ^ 

G(A*;e-ry) = EE ^ ^ 1 - C W > 0. (3.2) 

Observing that ry i—A*G(A*;^ • ry) only varies by the constant — ^ from the Dirac distribution 
motivates the following theorems which express a sufficiently smooth function by its integral mean 
value and a correction term involving Green’s function. For more details, the reader is again 
referred to m and the references therein. 
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Theorem 3.1 (Fundamental Theorem for A*). Let F be of class Then 

FiO = ^ / Fir^)duj{r,) + f G(A*; ^ ■ r,)A;FM dco^, ^ G O. 
Theorem 3.2 (Fundamental Theorem for V* and L*). Let F he of class C'^^^(n). Then 
FiO = ^I^Fir,)du;{ri)-I^A;GiA*-i-7j)-A;F{r,)dio{7j), ^ G 

where A* denotes one of the operators V* or L*. 


These theorems directly yield simple integral representations for solutions to the spherical differ¬ 
ential equations with respect to V*, L*, and A*. 

Next, we present a spatial regularization of G'(A*; •) around its singularity. This is a crucial step 
for the later definition of the scaling and wavelet kernels of the multiscale representation. 


Definition 3.3 (Regularized Green’s Function). Let Rp, p > 0, be of class C^"'^([— 1,1]), n G N 
fixed, satisfying 


lim 

P^O-l- 




dt = 0, 


fc = 0,1, 


and 



Then the function 



G(A*;t) 


t—l — p 


fc = 0,1,..., n. 


G^(A*;e-p) 


G(A*;^p), l-^.p>p, 


is called regularized Green’s function (of order n). RP is called the regularization function. 


A typical choice for Rp is the Taylor series of G(A*; •) centered at 1 — p and truncated at some 
power n. An exemplary plot for different scaling parameters p can be found in Figure Similar 
regularizations, but only for Taylor polynomials up to degree 2, have been used in other areas of 
geosciences, e.g., in and m- To be able to state a multiscale decomposition, it has to 

be guaranteed that convolutions with the regularized kernels converge to convolutions with the 
original kernels. The proofs are based on the fact that p G(A*;^ • p), p V|G(A*;^ • p) and 
p I—>■ L^G(A*;^ • p) are integrable on the sphere SI, uniformly with respect to ^ G SI, and can be 
found in [2] and [TT] . 

Lemma 3.4. Let Gp{A*-, •) be of class G^^)([—1,1]) and F of class G*^°^(S1). Then we have 


lim sup 
Jen 


G'’(A*;C-p)F(p)dcc(p)- [ G{A*-^-p)F{p)doj{rj) 
! Jq 


= 0 . 


Lemma 3.5. Let F be of class andGP{A*]-) of class 1,1]). Then 


lim sup 


A|G'’(A*;C-p)F’(p)dcc(p)-A| / G(A*; ^ • p)F(p)du;(p) 
! Jn 

where A* denotes one of the operators V* or L*. 


= 0 , 


Relations for higher order derivatives are simple consequences of the above lemmas by use of 
well-known surface versions of Green’s formulas that shift the differentiation from the convolution 
kernel to the convolved function F. Thus, they also require a higher smoothness of F. 
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I -original . rho=l/2 - rho=l/8 - rho=l/32| 


Figure 1: Plot of a twice continuously differentiable regularization Q G'^(A*; cos(0)) (left) and 
a twice continuously differentiable regularization Q ^ Sp{cos{6)) (right) at different scales p. 

Corollary 3.6. Let {/S.*] ■) be of class and F of class Then 

A*^GP{A*;i-p)Fir,)du;{p)-A*^ [ G{A*-f-p)F{p)du^{p) 

1 Jo. 

Corollary 3.7. LetGP{A*]-) be o/class C(^^([—l, 1]) and f of class . Then 

{{Al)^^{A;)^GP{A*;^-p))f{rj)dco{p) 


lim sup 


= 0 . 


lim sup 


-(At)^ {A*)^GiA*-,f- p) ■ f{v)dojiv) 


= 0 , 


where and A 2 denote one of the operators V* or L* (® denotes the tensor product x^y = xy"^, 
for x,y G 


An adequate choice of admits an explicit statement on the convergence rate. More precisely, 
if \RP{t)\ dt = 0{p) and |^i?^(t)| dt = 0(1), one can find 


f GPiA*;^-p)FiTj)du;{p)- [ G{A*-^-y)F{p)duj{p) 
in Jn 

[ AlGPiA*;^-p)Fip)du;{p)-Al [ G{A*-i-y)F{p)duj(,p) 
In Jn 

for F of class C*^°^(r2). If F is of class C^^^(f2), it even holds 


/ AlG^iA*-f-p)Fip)dujiy)-Al / GiA*;^ ■ p)Fip)dioip) 

In Jn 


= 0{pln{p)), 

= o{p-^), 

= 0{p\n{p)). 


The conditions on the regularization function are satisfied, e.g., by the choice of R^ as the truncated 
Taylor series of G(A*; •). 


3.2 Single Layer Kernel 

By the singel layer kernel we denote the convolution kernel of the integral operator with D 


formally given as in (2.8). Observing that A*Yn^k = —'n{n+^)Yn,k, the fractional pseudodifferential 
operator D, mapping the Sobolev space i7s(n) into can be defined via 


DF =(-A 


(-^‘+9 


_ 00 2n+l 

n—O k—1 




(3.3) 
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for of class where {■,-)L^{n) denotes the inner product (-F, G')l 2 (q) = J^F{r])G{ri)duj{r]). 

Its inverse D~^, mapping Hs_i(n) into Hs(fl), is correspondingly given by 


oo 2n+l 


D-^F= (- 


(^+ 4 ) ^ , 1 i^’^n,k)L^{a)Yn,k, 


(3.4) 


n—0 k—1 


for F of class Hs-i{^). From the addition theorem and the power series in Lemma 2.2 it is easy 
to derive the integral representation 


OO 2n+l ^ 1 /* 1 

—Fii^^Yn,k)L^{n)Yn,kiO = Tr-7^ / . F{r))duj{r]), ^ G fl. 


n=0 k=l 2 


The function 

^ ^ ^ , l-C-?7>0, (3.5) 

\/2 Vl - ■ ?7 

is from now on called the single layer kernel^ and denotes the starting point for our further con¬ 
siderations. The integral operator D~^ is called the single layer operator. For a more general and 
detailed overview on spherical pseudodifferential operators and the definition of Sobolev spaces, 
the reader is referred to, e.g., m and [26) . Since we are dealing with continuously differentiable 
functions in the remainder of this paper, it should be remarked that D~^ actually maps C^^^(n) 
into CW(ff), fc G No. 

In analogy to Green’s function, one can define a spatial regularization of the single layer kernel. 

Definition 3.8. Let p > 0 and a non-negative function of class C'^"^([—1,1]), n G N fixed, 
satisfying 


lim p'” f 
p^o+ 



dt = 0, 


A: = 0,1, 


and 



Then the function 



fc = 0,1,..., n. 


■ V) 


-R^(C ■»?), 1 - C • ?7 < P. 


is called regularized single layer kernel (of order n). 


The regularizing function Rp is generally chosen as the Taylor series of S centered at 1 — p and 
truncated at degree n. An exemplary plot for different scaling parameters p can be found in 
Figure The special cases of a linear or quadratic regularization have been applied to multiscale 
methods in physical geodesy, e.g., in |T0], [14]. In the Euclidean space the kernel S can be 
related to the fundamental solution of the Laplace operator A. A different kind of regularization 
for that kernel is treated in [T]. In our setting, we obtain the following limit relation in the same 
manner as for the Green function case in the previous subsection. 

Lemma 3.9. Let be of class 1,1]) and F of class Then we have 


lim sup [ ■ r])F{r])duj{r]) 

Jq 


S{^ ■ r])F{r])duj{r]) 


= 0 . 
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For the relations involving the surface gradient and the surface curl gradient, one has to observe 
that 77 !-)■ V|S'(^ • T]) and 77 i-)- • 77 ) are not integrable on the sphere fl. However, if F is of 

class and tj G denotes the rotation matrix with = ( 0 , 0 , 1 )^, we obtain 


S{^-r])F{'q)duj{r]) = V| / S{^ ■ tJr])F{tJr])du}{r]) 

! Jn 


= / 5 ( 773 )i^(tf77)^07(77) 

Jn 

= f ^(773)V|F(t[77)dc^(77), 

Jn 

where 77 = ( 771 , 772 , 773 )^ S H. Furthermore, regularizing the single layer kernel yields 


• v)Fiv)du;{v) = V| [ ■ 77 )^( 77 )^ 07 ( 77 ) 

1 JQ 


= VJ 


SP{r]3)F{tJr])du}{r]) 


= / 5'’(773)V^F(tf77)da;(77), ^ 

Jn 

so that the previous lemma implies the desired relations. 

Lemma 3.10. Let F be of class (7*'^^(H) and of class C^^^([—1,1]). Then we have 


lim sup 

p->-0+ 


■ 77 )^( 77 )^ 07 ( 77 ) - A| / Sif ■ 77 ) 1 ^( 77 )^ 07 ( 77 ) 


= 0 , 


where A* denotes one of the operators V* or L*. 


Relations for higher order differential operators follow analogously. Of more interest to us are 
combinations of the single layer operator with Green’s function for the Beltrami operator. It 
holds, e.g., that 


lim sup 




; e • 77 )) F{Tj)du;{p) - ^ G(A*;^ il)F{p)du^{r^) 


= 0 . 


(3.6) 


However, it is difficult to explicitly calculate ^G^(A*;^ • 77 ), as it would be required for our 
later applications. I?^^G(A*;^ • 77 ), on the other hand, can be calculated, and a regularization 
afterwards yields a similar limit relation. 


Lemma 3.11. For G LI we have, 


i? 7 lG(A*;^ 77 ) = —In (l + C- 77 ) -- 


27r 


2 1 - 25(e • 77 ) 


1 

27r 


Proof Lemma 2.1 and the pseudodifferential representation (3.4) imply 

-Pnif ■ v) 




1 277+1 


1 


—(n+i 47r —77(77 + 1) 

— i ^ 

00 1 °° 1 

"d')~ ■ d) 

ZTT ' 77+1 ZTT ^' 77 

n—1 n—1 

277 I 


Tf,) ) 


1 

27r 


ln((l + C-77)(^2 l-25(^77))) 


1 


which is well-defined for every 77 S H. 


^(l + ln( 2 )) 


□ 
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The above derived representation implies that (f, r;) i-t • ry) is zonal and of class 

X n). Some lengthy but basic computations yield 

vpyG(A-;«-„) = 1 (1 -SK.(,- «■ ,)f). C, e a 

A further application of the surface gradient causes a singularity of type Therefore, 

we do the following regularization for p > 0, 


4.K.1) = S'-K ■ 1) - 2 + -)■ 

For this kernel we can calculate 




(3.7) 


VJ®,»,.(,,{) = ^ ^ ^^ ) v| a a - K. ,)„) 

(-(*')'«■ V) + )!))» ) ® <«- «■ 

where denotes the one-dimensional derivative of . The analogous procedure works for the 
surface curl gradient, and we have for p > 0 that 


= i{l- S'« ■ 1 ) - 2 + 4S>k - ,) 




(3.8) 


and 




+T(-(s.)'({.,) + 


4(g^)'(g-p) 

(2 + 4,SP(C-p))2 




Thus, relation (3.6) can be formulated in the following numerically more advantageous way. 
Lemma 3.12. Let F be of class {Lt) and S'’ of class 1,1]). Then we get with s^,{-,-) 

{ . /-to nnrt 9 O ) /-im/] Q Q 1 nr-/-) onrt/->/-!■{■/i nt nln t ■f-i'i/ft- 


and sf,{-,-) as in (3.1) and (3.8), respectively, that 


lim sup 


SA-{^,v)Fiv)du}{Tj)-( G{A*-,f-r])F{r])du}{r]) 
1 Jn 


= 0 , 


where A* denotes one of the operators V* or L*. 


The relation we are actually aiming at, and which we require in later applications, is the following 
tensorial one. 

Lemma 3.13. Let f be of class c^^^(r2) and of class G^^^([—1,1]). Then we get with Sy.(-,-) 
and sf,{-,-) as in (3.1) and (3.8), respectively, that 


lim sup 
p-*-o+ {en 


(A| 0 s^. (p, f)) flji) duj{T]) ^G(A*; ^ • p)) • /(p) duj{r]) 


= 0 , 


where A* is one of the operators V* or L*. 
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Proof. Since |Vj 0 (C — (C ’ v)v)\ and |/(? 7 )| are uniformly bounded with respect to ^,r] G by 
some constant M > 0, we get the following estimate for ^ S and p > 0, 


(V| 0 (p, 0) fiv) Mv) - ^ (^€ ® e • v)) fiv) dco{v) 


(3.9) 


< 




■ v) + . , .L -v)- ^ 


2 + 4S{f-v) 

X IVf (g) (C - (^ • v)v)\ \f{v)\duj{r]) 


2 + ASP{^-r,) 


f SIC n)> (S-Vfr r) I 

/.=» (2 + 4SK-,))2 < ’'>+(2 + 4S<>({.,))= 




< 


X l(»? - (^ • ? 7 )^) v)v)\ \fir])\dujir]) 

1 . 1 


i-€ •'p<p 


*5(C • V) 


2 + 45(e • v) 


- ■ v) - 


2 + 4SP{^-r)) 


+M f S(f - _ (S^v. „) + GMiiA 

' (2 + 4S(f:„))2 < ’''+(2 + 4SP({.,) 


duj{r]) 

1 

vW 


Observing 


x\v-{^-v)i\\i-ii-v)r]\d(^iv)- 
\v - - {^ ■ v)v\ = ^ 


2S{^ ■ 77)2 




the integrability of 77 S{£, • 77 ) on the sphere O, and the properties for Sp from Definition 3.8 


we see that the integrals above vanish as p tends to zero. Due to the zonality of the kernels, 
this convergence is uniform with respect to ^ G 0. Furthermore, the convergence of (3.9) to zero 
additionally yields 

[ (v^OV;D-iG(A*;e- 77 ))/( 77 )dcc( 77 ) = V^ f (v;D- 1 G(A*; ^ ' 77 )) • 7 ( 77 ) ^ 02 ( 77 ), 

and therefore, the desired statement. The assertion for the surface curl gradient follows analo¬ 
gously. □ 

In analogy to the Green’s function case, an adequate choice of RP admits an explicit statement 
on the convergence rate. More precisely, if \RP{t)\ dt = 0{p) and dt = G(l), 

one obtains 


• v)Fip)du;{p) - Af / ■ 77 )^( 77 )^ 02 ( 77 ) 


= o{p-^). 


For the combination of Green’s function and the single layer kernel the same convergence rate 
holds true. 


(A|®s^.(77,^))/(77)dw(77)-A^ f (^A*Dj^G(A*;^- 77)) •/(77)da;(77) 

J Q 


= oip-^), 


for F of class G*-^^(n) and / of class c*^^^(D). The conditions on the regularization function are 
satisfied, e.g., by the choice of RP as the truncated Taylor series of S. 


4 Spherical Decompositions 

We introduce two decompositions relating to the two sets of vector spherical harmonics from 
Section 2 The first one, relating to y^\ and the operators respectively, is the well-known 
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spherical Helmholtz decomposition. It decomposes a vector field into its radial part and two 
tangential parts. The second one, relating to and the operators 0 *-*^ respectively, is crucial 
for the separation with respect to the sources and is presented in some detail in this section. A 
more general overview on similar spherical decompositions can be found, e.g., in |15j . 

Theorem 4.1 (Helmholtz decomposition). Let f be of class Then there exist uniquely 

defined scalar fields Fi of class andF 2 ,F^ of class satisfying 

J^Ffir])duj{r]) = 0, z = 2, 3, 


such that 


/(C) = ofVi(C) + of A2(C) + ofF3(C), Cea 

A proof of this decomposition can be found, e.g., in [3]. Using Green’s function for the Beltrami 
operator and Theorem |3.2| yields the following representations for the Helmholtz scalars. 


m) 




(4.1) 

m) 

= - / (v;g(a*;c-^)) 

Jo. 


C s u, 

(4.2) 

m) 

= - [ {L;G{A*;C-y))- 
JQ 

fiv)duj{r]), 

c G II. 

(4.3) 


If Fi additionally has vanishing integral mean value, i.e., J^Fi{r])duj{r]) = 0 (as is the case for 
functions satisfying the pre-Maxwell equations), there exists a function U of class with 

A*U = Fi, such that Theorem 3.1 implies 


Fi(C) = A| / G(A*;C-^h-/(^)dw(^), 


(4.4) 


While the orthogonality is the main property of the Helmholtz decomposition, a representation 
with respect to the operators is of special interest in geomagnetic modeling. In order to obtain 


a representation of the corresponding scalars, we rewrite (2.5)-(2.7) as 



2 2 

(4.5) 



(4.6) 

o(3) 

= 5(3). 

(4.7) 


This gives us the necessary representation to prove the following decomposition theorem. 

Theorem 4.2. Let f be of class Then there exist uniquely defined scalar fields Fi,F 2 of 

class G^^^(H) and F^, of class satisfying 

^ [ h{r])duj{ri) = 0 , 

j r2 

Y [ h{r]) - h{jg)duj{r]) = 0 , 

J Q 

such that 

/(c) = df^F^{o+5fm)+ofm), Cell. 
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The scalars Fi, F 2 and F 3 can be represented by 


Fi 

= -D-^Fi + 174 - 1+2 - -F 2 , 

2 4 2 ’ 

(4.8) 

F 2 

= Id-^f, + Id-^F2 + If2 , 

(4.9) 

F 3 

= F 3 , 

(4.10) 


with Fi , F 2 and F 3 being the uniquely determined functions of the Helmholtz decomposition in 


Theorem 4-1 


Proof Applying the Helmholtz decomposition to / and using (4.5)-(4.7), we get on 

= ^ 5 ^^) D -^ Fi + ^5(2)74-17^1 + ^ 5 ( 1 ) Q74-1 - 1^ F2 + ^5(2) Q74-1 + 1^ A2 + 5(3)7^3 

= 5 ( 1 ) + \d-^F 2 - ^^ 2 ^ + 5 ( 2 ) (^Id-^Fi + \d-^F 2 + ^F 2 ^ + 5 ( 3 ) 7 ^ 3 . 


This implies a decomposition as stated in the theorem. Due to the uniqueness of the Helmholtz 
representation, it follows directly that F 3 is defined uniquely when having a vanishing integral 
mean value. For the uniqueness of Fi and F 2 it is sufficient to show that /(^) = 0, ^ € H, only 
has the trivial decomposition with respect to the operators o^*), i = 1,2,3. If 


+ 51 ^ 1 ^ 2(0 + dfF.iO =0, e e H 


A2) 


(3), 


we get from (2.5)-(2.7) that 


di) 


(( 74 ^ + i) - i) Mo) + {f 2{0 - F,{o) + ofF^iiO = 0 , e e 


The uniqueness of the Helmholtz decomposition then implies 

F2{0-M0 = 0, 
(745 + 1) A (0 + (^«-2)^2(e) = 0, 


if 4 F In Fiiv) - Mv)doj{r]) = 0, which gives us 


DMii 0 = 0 


Thus, F’i(C) = 0, ^ G H, since 74 is injective, and it follows F 2 (C) = 0, ^ G H, so that uniqueness is 
given for this decomposition. □ 


The theorem above yields, by use of (4.1 )-(4.3), a representation of the scalars Fi. Of importance 
in the later application, however, are the vectorial quantities o^'^'^Fi. Thus, we first calculate from 
(4.8) that 


O^^MO = • fiO) + ■ /(C)) - ^ v;g(a*;C • v) ■ fiOdojiO 

+ Ud^ j V;G(A*;C-^)-/(^)c?cc(r,) (4.11) 

J fi 

-^V|74/i(C-/(C)) + ^V|74/i [ V;G(A*;C-r?)-/(7?)dcc(77) 

J fi 
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The expression in the second row, involving the operator D, is unfortunate since we have no 
explicit representation for the corresponding regularized convolution kernel. Observing 

this can be circumvented by rewriting 


JQ. 


-1 


= 4^^ 


V;G(A*;e-ry)./(,7)dcc(r;) + i?^-iA| [ GiA*;^ ■ ■ ftan{v)doj{v) 

! Jn 


j V;G(A*;C • v) ■ fiv)Mv) + ' ftaniO, ? e O, 

J Q 


3.1 


has been used in the last step, assuming / to be of class Application 


where Theorem 
of the above to (|4.11|) provides an easier calculable expression 


• fiO) + ■ m) + ■ ftUO - ■ fiO) 

+ i V| V;Oj-iG(A*; e • v) ■ fiv)duj{n) - ^ V;G(A*; ^ r,) • fiv)dio{v). 

Since the occurring differential operators and the integration cannot be interchanged without 
restriction, we need to switch to the regularized versions of the single layer kernel and the Green 
function for the Beltrami operator. Then it is valid to set 

• /(O) + ^ fiv)duj{v) - V;5'^(C • v) ■ f{v)duj{r]) 

~4~ / fi'n)di^iv) + \ f s^4v,0) fiv)duj{r]) (4.12) 

»/ i/ O 


[ (V^OV;G^(A*;e-r;))/(,7)da;(r?), ^ 

Jn 


for p > 0. If G^(A*; •) is of class G^^^([—1,1]) and of class G^^)([—1,1]), the considerations in 
Section imply 


lim sup 

P-J-0+ JgQ 


/W(e)- 6 WA (0 


= 0 . 


(4.13) 


Analogous computations for F 2 yield a regularization 

= |C(C •/(?))-/" SP{^-r])iT f{r])duj{r]) +f V*SP{^-T])-f{T])duj{r]) 

+ ^ / VjS'^(C-7?)?7-/(?7)da;(77)-^ / {^1 ® s^,{r],^)) f{r])duj{r]) (4.14) 

1 / r 2 r 2 

^ (v| 0 V;G^(A*;e • v)) f{ri)d 0 Jiv), ? € n, 


lim sup 

p->- 0 +^gQ 


/f (0-5(2) F 2(0 


= 0 . 


(4.15) 


Finally, the determination of corresponds to the calculation of the toroidal part of /. From 

(4.10) and (4.3), we get 

5 ^^^hiO = -Li[ L;G{A*;C-ri)-f{rj)dojir,), ^ G n. 

Jn 
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Corollary |3.7| then implies for the regularized version 

/f (0 = - / {LI 0 r,)) f{rrj)dLu{r,), ^ G ft, 

Jn 

that for G^(A*; •) of class G(^^([—1,1]), 


lim sup 

p-i-0+ JgQ 


/f)(e)-5(3)F3(e) 


= 0 . 


(4.16) 


(4.17) 


Summarizing, we can state the following theorem. 

Theorem 4.3. Let / be of class c*-^^(n), with J^rj ■ f(r/)duj(r]) = 0. Furthermore, let the regu¬ 
larized Green function G^(A*;-) be of class G^^^([—1,1]) and the single layer kernel of class 
G(i)([-l,l]). Then 


fif) = 5(i)Fi(e) + 5(2)F2(e) + 5(3)F3(e), C G 


with 


lim sup 

p->-0+ JgQ 


^p\^,v)f{v)dujir]) - oW|;(C) 


= 0 , 


for i = 1, 2, 3. The convolution kernels are given by 

+ |V| 0 (r?, 0 - ^V|^^(e, 7?) 0 77 - 1 V| 0 \7;GP{A*;f ■ 77 ), C, ^ G fl, 

= ^®77QA*G^(A*;e-77)-^5'’(C-77)) +^^®V;.S^(C-r7) 

-iv| 0 ( 77 ,0 + ^V|^^(e, 77 ) 0 77 - 1 V| 0 V;G^(A*; e • 77 ), C, ^ G fl. 


^f(e,^) = -G|0T;G^(A*;e-77), C,^Gfl. 


Proof Since f^r/ ■ f(g)duj{r]) = Fi(ri)duj(ji) = 0, representation (4.4) implies 


^G(A*;C-77)77-/(r;), 

Substituting the Green function by its regularized counterpart, we obtain 
^^A^ ^G^(A*;^-77)77 •/(77)da;(?7) = ^ A^G'’(A*; f • 77)77 • 7(77) ^07(77) 

= ^ QaJG^(A*;C- 77)^^/(?7)da;(77), ^ e fl. 


Analogously, the remaining integral expressions in (4.12) can be written in terms of convolutions 
with tensorial kernels if this is not already the case. This yields the kernel $^^^(-,-) for fl^\ 


Corollary 3.6 and (4.13) provide the desired limit relation for The same holds true for 

d^^^F2 and^^^3. □ 


5 Multiscale Representation for the Separation of Sources 


The kernels from Theorem 4.3 are the main ingredient to the upcoming multiscale representation. 
They actually denote the so-called scaling kernels, while the differences for different parameters p 
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denote the corresponding wavelet kernels. But before we go into detail, we briefly want to motivate 
why the decomposition with respect to the operators 5^*^ can be called a separation with respect 
to the sources. 

From now on, we denote the magnetic field by b of class the corresponding source current 

density by j of class and by fiQ we mean the vacuum permeability. Furthermore, we 

assume the pre-Maxwell equations 

\7^Ab{x) = fJ-oj^x), X G K^, 

Va: • b{x) = 0, X G K^, 

to be satisfied. As mentioned in the introduction, the Mie decomposition (see, e.g., [5] and 0 ) 
yields poloidal fields pb, Pj and toroidal fields qb, qj, such that 

b{x) = pb{x) + qb{x), X G 
j{x) = pj{x)+qj{x), xGK^. 

Making use of the law of Biot-Savart (see, e.g., [12]) and the fact that the poloidal magnetic field 
Pb is solely produced by tangential toroidal current densities qj , the poloidal magnetic field can be 
split up as follows, 


Pbix) =pI^*{R;x) +pI^ 

where 

W,Apl^\R;x) = 
y,ApT\R-x) = 


In other words, _p™*(i?;-) denotes the part of the magnetic field that is due to source currents 
in the interior of the satellite’s orbit 11^, and •) the part due to source currents in the 

exterior. A more detaile d des cription can be found, e.g., in m and [22]. Since •) is still 

divergence-free, equation (5.1) implies that a harmonic potential C/“‘(i?; •) : —)■ K exists, such 


x), X € 

E \ flfi. 


Poqjix), 

0, 

X G 

X G 

(5.1) 

0, 

X G 

(5.2) 

Poqj{x), 

X G 


that p™*(i?;x) = Va;[/*"‘(i?;x), for x G The potential [/“* can be expanded w ith re spect 


to the outer harmonics and the application of the gradient in combination with (2.111 then 


implies that p™^{R- •) can be expanded in with respect to yn\- Analogously, pl^*{R; •) relates 

Vn k- remaining toroidal part qb can be interpreted 


to an expansion in 11)^* with respect to 


as the part induced by poloidal source currents Pj crossing the sphere Vtn, and corresponds to the 

~(3) 

vector spherical harmonics To sum up, additionally observing that the poloidal and toroidal 
fields are continuous up to we find 


b{x) =pI^\R-x) 


-pI^\R;x) + qbix), X G flfl, 


(5.3) 


with yll\, i = 1, 2,3, n G Ng., fc = 1, ..., 2n-|-l, being the appropriate basis system for this split-up. 


^n,k 


Finally, the definition of the vect or sp herical harmonics in (2.9) implies that Theorem 4.2 yields 
the exact same decomposition as (5.3). More precisely, the o'^part denotes the contribution due 
to sources in the o^^^-part the contribution due to sources in and the o^^^-part the 

contribution due to source currents crossing the sphere fl/j. 


5.1 Multiscale Representation 


Now, we turn to the actual multiscale representation. We discretize the regularized kernels from 
Theorem 4.3 by choosing parameters p = for J G Ng. The scaling kernels (of scale J) are 
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Figure 2: Absolute value of the scaling kernels (left) and wavelet kernels (right) for the interior 
poloidal contribution at scales J = 1,4,8 (note that the color scaling is logarithmic). 

then defined by 


II 


e fl, 



e fl. 



^,77 G fl. 


These kernels still have global support. The announced locally supported wavelets are obtained 
by taking the difference of two such scaling kernels. A wavelet kernel (of scale J) denotes one of 
the following kernels 

Due to the regularization of the Green function and the single layer kernel, these wavelets clearly 
have local support in a spherical cap of radius 2“"'. More precisely, we find that supp(’®'j(^, •)) C 
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{t] G ^\1 — ^ ■ rj < 2 for i G {int, ext, q\ and ^ G An illustration of the kernels is given in 
Figure 

Each scaling kernel generates a scaling transform. These are given by 

P}6(i) = [ ^:j{^,fj)b{Rfj)dojirj), x = R4GnR. (5.4) 

Jn 

The corresponding wavelet transforms read 

R:jb{x) = f ^){^,r,)b{Rr,)duj{r^), x = R^ G Hr. (5.5) 

Jn 

The idea of the multiscale approach is to resolve the modeled quantities at different spatial resolu¬ 
tions. Therefore, the scaling kernels are only used to provide a trend approximation of the coarse 
features at some small initial scale Jq. The spatially stronger localized features are subsequently 
added by the wavelet transforms. This is reflected in the following relations, 


.7-1 

P}b{x) = P}^b{x) + Y, x = R^G Hr, (5.6) 

j=Jo 

for i G {int, ext, q}. Different from the multiresolution constructed, e.g., in [52] and [23], the scale 
spaces V} = {P}b\b G c^^^(Dfl;)} in our approach are not necessarily nested in the sense V} C 
The advantage here is the local support of the wavelet kernels, which implies that the evaluation 
of the wavelet transforms R{jb{x), for x = R^ G Hr, only requires data in a spherical cap around 
^ G H with scale-dependent spherical radius 2~^. Thus, regions with a higher data density can be 
resolved up to higher scales, i.e., up to a higher spatial resolution, without suffering errors from 
the lower data densities in surrounding areas. The general concept is illustrated by the following 
tree algorithm 



The maximal scale J = Jmax at which Pjb{x) can be evaluated is determined by the amount 
of data points in the vicinity of x = R^ G Hr. Sufficiently many data points in the support of 
are required to guarantee a numerical meaningful evaluation of the integral in the 
wavelet transform i?j _ib{x). 


To conclude this subsection, we summarize the results in the following theorem, which is mainly a 
reformulation of Theorem 4.3 in terms of the above described multiscale setting for the separation 
of the magnetic field with respect to its sources. 


Theorem 5.1. Let b be of class c*^^)(IR^), and j of class c^^^(M^), satisfying the pre-Maxwell 
equations 


V2;A5(x) = p-o3{x), xG 

Vx ■ b{x) = 0, X G 


If PT\ P' 


J ’ ^ J’ 


. RT, p'!j are defined as in |f5.5l), then 

b{x) = p“‘(i?;x)-bp6'^‘(i?;x)-f gh(x), x G Hr, 
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for a fixed R > 0, with 


pl-\R;x) 

oo 

= P}"‘6(cr) + 5] P“‘6(x), xenn, 


3=Jo 

pT\R\x) 

OO 

= Pjf 6(x) + ^ Pf‘6(x), xG^e 


j=Jo 

qb{x) 

= Plbix)+f2R'jb{x), xG^n. 


3=Jo 


5.2 Crustal Field Modeling from CHAMP Data 

In this subsection, we apply the above derived multiscale approach to a set of CHAMP satellite 
measurements. The used data set is similar to the one used in [22] and |23| and has been collected 
between June 2001 and December 2001. It has been pre-processed at the GFZ Potsdam by Stefan 
Mans to fit the purpose of crustal field modeling (see, e.g., HI] for a detailed description). Due 
to the almost spherical orbit of the CHAMP satellite, we can assume all data to be given on a 
sphere of radius Re + 450km, where Re = 6371.2km denotes the mean Earth radius. For the 
discretization of the integrals appearing in the multiscale representation, we use the integration 
rule described in |6], which requires an equiangular data grid (and reflects the data situation of 
satellite measurements). In this example, we use a grid with 180 grid points in latitudinal as 
well as longitudinal direction. Centered around each grid point, we select a spherical rectangle 
with a diameter of 2.5° in latitude and longitude and average all measurements in the cell using a 
M-estimation with Huber’s weight function (see, e.g., HHl). The resulting input data set is shown 
in Figure 

The results obtained from the multiscale representation are illustrated in Figures |4]|^ For the sake 
of brevity, we only indicate the radial and the south-north component of the magnetic field, and in 
Figure]^ only the radial component. Figure]^ also illustrates best the different spatial resolutions 
of the multiscale representation. The initial trend approximation at scale Jq = 2 only resolves 
very coarse features, while the subsequent wavelet transforms resolve more and more localized 
features, such that the scales J = 5,6, 7 mainly focus on the strongest crustal field anomalies over 
Central Africa and Eastern Europe, as well as North America and Australia. In oceanic regions, 
there is hardly any contribution at these scales, indicating that there the crustal field is of a rather 
coarse nature (i.e., of large wavelength when arguing in frequency domain). Eurthermore, one 
finds that the structure of the wavelet contributions hardly changes for scales higher than J = 5. 
This might be an indicator of the general spatial extend of the anomalies of the crustal field signal 
at satellite altitude (when comparing the resolved features with the size of the support of the 
wavelet kernels). Eigure shows the final approximation P}"* b of the internal magnetic field 


I 

t 




Figure 3: The radial component (left) and the south-north component (right) of the input magnetic 
field (in nT), averaged to a 180 x 180 equiangular grid. 
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Figure 4: Approximation of the internal poloidal magnetic field (in nT), at initial scale J = 2 
(top left) and at scale J = 8 (bottom right); only the radial component is shown. The remaining 
figures show the intermediate wavelet contributions from scale J = 2 to scale J = 7 (top right to 
bottom left). 
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Figure 5: Approximation of the internal poloidal magnetic field (in nT), at scale J = 9; radial 
component (left) and south-north component (right). 




Figure 6: Difference between the input data set and the approximation of the internal poloidal 
magnetic field (in nT), at scale J = 9; radial component (left) and south-north component (right). 


contributions at the highest scale Jmax = 9. The difference to the input data set is indicated in 
Figure It actually illustrates the performance of the separation with respect to the sources. 
One can recognize strong polar fields that are clearly not due to the Earth’s crustal field and are 
probably induced by polar ionospheric current systems. Furthermore, one finds bands of positive 
and negative field strength oriented parallel to the dipole equator. This is a typical signature of 
magnetospheric ring currents. 

Thus, the multiscale approach of this paper can be used to improve pre-processed crustal magnetic 
field data. The decomposition with respect to is actually able to filter out contributions 
originating outside the satellite’s orbit or at satellite altitude, while the wavelet transformations 
give a clearer impression of the local features of the crustal field. 
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